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1.  Introduction 


It  has  been  shown/fSJfl^that  modern  DSP  (Digital  Signal  Processing)  and 
Pattern  Recognitiorntgcnniques  are  very  effective  tools  for  Ultrasonic  NDE 
(  Non-Destructive  Evaluations  )  of  materials.  In  our  recent  researches  [1], 
such  techniques  as  High  Resolution  Spectral  Analysis,  Deconvolution, 
K-mean  Clustering,  Non-linear  Mapping  and  Multipe  Fisher’s  Discriminant, 
have  been  applied  to  the  classification  of  hidden  defect  geometries  with 
high  accuracy. 


In  view  of  this  uprising  needs  for  DSP  and  Pattern  Recognition  in 
Ultrasonic  NDE,  a  software  package  that  effectively  combines  the  two  is  in 
order.  Interactive  Ultrasonic  NDE  ( lUNDE  )  developed  by  Information 
Research  Laboratory,  Inc.  is  such  a  package  created  especially  to  meet 
the  needs  of  researches  and  practising  engineers  in  this  field.  It  also 
includes  a  variety  of  feature  extraction  functions  that  allow  the  user  to 
choose  features  extracted  from  the  time  series  of  his  or  her  interest  for 
classificationrThe  po^er  of  this  package  is  further  enhanced  by  its 
versatile  graphics  dismay  capabilities,  which  can  produce  scatter  plots, 
2-dimensional  and  3-dimensional  graphics  on  a  monitor  and  on  a  plotter  or 


printer. 


In  fact,  the  DSP  and  Pattern  Recognition  libraries  in  this  package  is 
so  rich  that  it  can  be  used  either  as  a  signal  processing  or  a  pattern 
recognition  package  or  both  without  lim  iting  to  NDE  applications. 


2.  Software  and  Hardware  Requirements 


Software  requirement : 

1 .  MS  DOS  (R)  3.0  or  above. 

2.  The  NDE  package  is  developed  under  the  Microsoft  (R)  C  5.0 
environment  using  the  large  memory  model.  It  is  strongly 
suggested  that  the  same  compiler  (5.0  or  above)  be  used  in  case 
the  user  needs  to  link  this  package  with  his  own  applications. 


Hardware  requirement ; 

1 .  IBM's  PC/XT  or  PC/AT  or  compatibles. 

2.  Graphics  adaptor  (HGA,  CGA  and  EGA)  installed. 

3.  At  least  280  Kbytes  of  free  on-board  memories. 
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3.  Setting  up  the  package 


Before  starting  to  run  this  package,  the  user  should  first  set  up  his 
programming  environment  for  graphics  display  by  running  the  EQUIP.exe 
program.  This  program  will  create  the  file  CONFIG.gpc  which  specifies 
the  type  of  printer/plotter,  monitor  and  graphics  adaptor  that  the  user  is 
using. 


NOTE :  If  you  have  a  harddisk,  the  most  convenient  way  is  to  copy 
all  the  files  into  the  harddisk  and  set  the  harddisk  drive  as  the 
defauit  drive.  Hnyway,  the  files  CONFIG.gpc,  CONFlG.peR; 
SIMPLEK.mt  and  COMPLEK.fnt  must  be  in  the  default  drive. 
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4.  Function  Descriptions  and  Operation  Guide 


Besides  the  EQUIP.exe  program  which  is  used  for  package  set-up,  there 
are  3  executable  files:  DEMO.exe,  NDE.exe  and  WIG.exe. 


4.1  DEMO.exe 

DEMO.exe  is  a  demonstration  of  some  of  the  functions  offered  by  this 
package.  It  requires  to  read  3  data  files:  SAMPLE,  SA  and  SB.  It  takes 
about  2  minutes'  time  to  read  all  these  data  on  an  IBM-PC/XT,  so  "  Be 
Patient ! " 


4.2  NDE.exe 

NDE.exe  can  perform  1 5  kinds  of  digital  signal  processing  or  pattern 
recognition  operations.  The  available  options  are  tabulated  in  the 
MAIN_MENU. 

MAIN^MENU 

bsp  Power  Spectrum  by  Burg's  Technique 

near  Nearest  Neighbour  Classification 

c  r  Cross-Correlation  between  two  waveforms 

dis  Graphics  Display 

f  e  Feature(s)  Computation 

fsp  Magnitude  Spectrum  by  FFT 

hb  Find  Analytic  Part  of  the  input  signal  by  FFT 

nimp  Nonlinear  Mapping 

spi  Spiking  Filter  Deconvolution  (no  IBM-PC  version) 

wv  Wavelet  Transformation 

win  Wiener  Filtering  Deconvolution  using  FFT 

wind  Hann’s  frequency  windowing  (a  low-pass  filter) 

fish  Multiple  Fisher's  Discriminant 

foly  Foley-Sammon  Transformation 

kmean  K-mean  Clustering 
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Once  entered  the  main_menu,  the  user  is  prompted  to  provide  his 
choice.  The  user  can  choose  the  desired  operation  by  entering  the 
corresponding  mnemonics  in  small  letters  (e.g.  bsp  <return>). 

After  entering  an  option,  the  user  should  be  ready  to  answer  a  few 
questions  to  provide  sufficient  informations  that  are  necessary  for  that 
particular  operation  (e.g.,  where  to  get  the  input  data;  how  long  is  the  input 
data:  where  to  store  the  results:  etc.).  In  the  following  few  sections,  we 
are  going  to  discuss  each  of  these  operations  in  some  details.  In  most  of 
the  cases,  we  assume  the  input  data  to  be  read  from  file(s)  and  the  results 
to  be  stored  in  file(s).  The  user  will  be  prompted  to  provide  the  necessary 
filenames  in  these  cases. 

In  the  following  subsections,  we  will  describe  each  of  these 
functions  in  some  detail.  Unless  stated  othenwise,  x(n)  and  XQ(n)  are  used 
to  denote  the  Test  Signal  and  the  Reference  Signal  respectively  while  X(k) 
and  XQ(k)  are  used  to  denote  their  Fourier  Transforms. 


4.2.1  bsp  (fig. 3) 

description 

This  function  computes  the  power  spectral  density  of  the  input 
signal  by  Burg's  technique  [2].  This  is  a  model-based  spectral 
estimation  technique.  The  input  time  series  is  modeled  as  an 
autoregressive(AR)  process  driven  by  white  noise,  such  that 

P 

X(k)  »  1  /  (  1  +  X  a[i]exp(-j27iik/N)  ) 

i»1 


where  P  is  known  as  the  AR  order  and  the  a[i]  's  are  known  as  the  AR 
coefficients.  There  are  no  general  rules  for  choosing  an  optimal 
value  of  P.  As  a  rule  of  thumb,  this  can  be  taken  to  be  about 
one-third  of  the  record  length  N  for  small  N. 

The  computation  is  divided  into  2  steps.  First  we  compute  the 
a[i]'s  using  Burg's  Technique  which  is  a  very  popular  one  among  some 
other  possible  methods.  Then  we  compute  |X(k)p  using  FFT. 
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informations  required 

1 .1  Filename  of  the  input  data  (e.g.  T 1  SaO.col). 

1 .2  Starting  point  and  ending  point  of  data  (e.g.  100,355). 

1 .3  Length  of  the  input  data  (e.g.  256). 

note:  maximum  is  512 

2.1  Desired  order  of  the  AR  model  (e.g.  20). 

2.2  Number  of  spectral  points  required  (e.g.  200). 

note:  This  is  a  resolution  issue.  The  user  should  be  aware  to 
choose  a  sufficierafy  large  number  such  that  no  significant 
structures  of  the  power  spectrum  will  be  missed.  The 
maximum  allowed  is  512. 

3.1  If  you  want  to  store  the  AR  coefficients,  you  need  to 
provide  a  filename  to  store  these  values.  The  number  of 
output  points  to  be  stored  in  this  example  should  be  20  ( 
equals  to  the  AR  order  used  ). 

4.1  An  output  filename  to  store  the  power  spectrum  (e.g. 
psd.dat). 

4.2  Number  of  output  points  to  be  stored  (100  for  this 
example).  This  is  because  the  input  signal  is  real  which 
implies  a  symmetrical  power  spectrum. 


4.2.2.  near 

description 

Suppose  we  have  c  classes  of  n  -vectors  in  an  n  -dimensional 
feature  space  and  that  each  class  contains  P  patterns. 

Let  t  a  the  test  vector  whose  class  identity  is 

of  interest  and  p  »  {pj,p2,...,pp)  be  one  of  the  known-class 
patterns  (or  training  pattterns).  The  Euclidean  distance 

is  computed  for  all  p  .  Then,  by  the  nearest  neighbor  rule,  the 
class  of  t  is  taken  as  the  same  as  that  of  the  p  with  the 
smallest  d. 
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in  form  a  tions  required 

1 .1  Number  of  features  in  each  pattern  (e.g.  3). 

1 .2  Total  number  of  training  patterns  (e.g.  12). 

1 .3  A  filename  for  each  feature  (i.e.,  3  filenames  needed  in  this 
example). 

1 .4  Starting  points  and  ending  points  to  read  these  features 
from  each  of  the  files. 

1 .5  Number  of  clusters  in  the  training-set  (e.g.  4). 

note:  For  this  example,  let  the  3  input  files  be  FileJ ,  File_2 
and  File_3.  Therefore,  File_} ,  File_2  and  File_3  should  contain 
respectively  the  feature(iys,feature(2ys  andfeature(3ys  of 
all  training-patterns.  In  other  words,  12  pieces  of  data  are  to 
he  read  from  each  of  these  files  in  this  example.  The  relative 
order  of  these  patterns' features  in  each  file  must  be  the  same 
and  that  we  have  assumed  that  there  are  equal  number  of 
training  patterns  in  each  cluster  (  3  training  patterns  per 
cluster  in  this  example  ). 


4.2.3.  cr  (fig.5) 

description 

This  function  computes  the  normalized  cross-correlation  of  2 
input  signals.  If  we  let  x(k)  be  the  first  signal  and  y(k)  be  the 
second  signal,  this  subroutine  will  compute: 

r^{m)/{  r„(0)r„(0)  )’'= 


where  '■xy(W  “ 

w 


informations  required 

1 .1  Filename  of  the  first  data  set. 

1 .2  Starting  point  and  ending  point  of  the  data. 

1 .3  Filename  of  the  second  data  set. 

1 .4  Starting  point  and  ending  point  of  the  data. 

1 .5  Length  of  each  data  set  (maximum  is  256). 

2.1  An  output  filename  to  store  the  result  of  cross-correlation. 
2.3  Number  of  output  points  need  to  be  stored.  This  would  be 
equal  to  2  times  minus  1  the  number  given  in  1 .5  above. 
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4.2.4.  dis 


! 

description 

Graphics  display. 


auailable  options 

(1)  scatter  plot,  or  (2)  waveform  plot. 


informations  required 
(i)  For  scatter  plots: 

1.1  Number  of  clusters  to  be  plotted. 

1 .2  Number  of  points  on  each  cluster. 

1 .3  Choose  a  symbol  to  represent  each  cluster  on  the  plot. 

note:  a  symbol  is  an  integer  from  0-7.  These  symbols  are 
respectively  shown  as  follows: 

□  QV  +  XOA* 


2.1  A  filename  to  read  the  data  to  be  plotted.  Each  piece  of  data 
should  contain  2  parts  -  the  (x,y)  coordinates. 

3.1  The  title  for  the  plot. 

3.2  The  label  on  the  x-axis. 

3.3  The  label  on  the  y-axis. 

(ii)  For  waveform  plots 

1 .1  Filename  of  the  input  waveform  data. 

1 .2  Starting  point  and  ending  point  of  the  data  to  be  plotted. 

2.1  The  title  for  the  plot. 

2.2  The  label  on  the  x-axis. 

2.3  The  label  on  the  y-axis. 

2.4  Number  of  data  points  to  be  plotted. 


At  the  end  of  each  plot,  the  user  can  return  to  the  main  menu  by 
hitting  the  <return>  key  or  he  can  dump  the  plot  onto  a  printer/plotter 
by  hitting  the  <m>  key  and  specifying  the  location  and  size  of  the  plot 
on  the  paper.  However,  the  user  may  suffer  from  a  shortage  of  system 
memory  that  forbids  him  to  dump  the  plot  onto  the  printer/plotter.  If 
this  is  the  case,  the  user  may  want  to  use  a  DOS  command  to  dump  the 
on-screen  plot  onto  a  printer  {  or  he  can  use  some  other  package  that 
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will  do  the  job  ). 


4.2.5.  fe 

description 

This  is  a  feature  extraction  subroutine.  The  several  options 
included  here  have  been  tested  to  be  very  effective  in  ultrasonic 
NDE  problems . 


auailable  options 

ar  Amplitude  Ratio. 

bk  Store  the  computed  features  and  return  to  the  main 

menu. 

bnw  Bandwidth  of  the  input  spectrum. 

kusk  Kurtosis  and  Skewness  of  the  input  time-series. 

fr  Frequency  Ratio. 

mcr  Maximum  of  Cross-correlation  with  reference, 
fkusk  Kurtosis  and  Skewness  of  the  input  spectrum, 
sft  Show  the  computed  features  in  tabulated  form. 

After  entering  an  operation  (in  small  letter),  again,  the  user  is 
prompted  to  answer  a  number  of  questions  to  supply  sufficient 
informations  that  are  necessary  for  that  particular  operation.  These 
operations  are  explained  below. 

(i)  ar 

description 

Amplitude  ratio  (ar)  is  defined  as: 

ar.{Z|x(k)|  }/{S|Xo(k)|} 
informations  required 

1 .1  Number  of  data  points  in  each  of  the  input  signals. 

1 .2  Filenames,  starting  points  and  ending  points  of  the 
referenced  signal  and  the  test  signal. 
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(ii)  bk 


description 

This  option  allows  the  user  to  store  the  computed  features  in  a 
file  and  returns  him  to  the  main  menu. 

informations  required 

If  the  user  wants  to  store  the  computed  features,  he  will  be 
prompted  to  provide  an  ouput  filename. 


(iii)  bnw 

description 

This  option  computes  the  bandwidth  of  the  test  signal  from  its 
magnitude  spectrum.  It  also  computes  the  frequency  of  peak 
magnitude  as  well  as  the  fractional  power  distribution  in  1 0 
equally  spaced  frequency  bands  (starting  from  d.c.  to  the  100th 
frequency  bin). 

note:  The  input  to  this  operation  should  be  the  magnitude 
spectrum. 


information  required 

1 .1  Number  of  spectral  data  points  available  in  the  input 
magnitude  spectrum. 

1 .2  Filename  of  the  input  spectral  data. 

(iv)  kusk 
description 

This  function  computes  the  kurtosis  and  skewness  of  the  input 
signal,  say  x(k),  where: 

Kurtosis  -  {  X^(k  -  \Lf  x(k)  }/ {  S^(k  -  p)2  x(k)  Y 
Skewness  =  {  X^(k  -  p  x(k)  }/ {  Z|^(k  -  |i)2  x(k)  Y'^ 

P-  k*x(k)/Z^x(k) 

1  0 


description 

The  definition  of  Frequency  Ratio  (fr)  is  the  same  as  that  of 
ar  in  (i)  above  except  that  |x(t)l  and  lxQ(t)|  are  replaced  by  lX(f)| 
and  |Xo(f)l  respectively. 

informetions  required 

Same  as  (i)  except  that  the  input  data  x(t)  and  XQ(t)  are 
replaced  by  |X(f)|  and  |XQ(f)l  respectively. 

(vi)  mcr 

description 

Let  a(k)  be  a  referenced  signal  and  x(k)  the  test  signal.  Define 
cross-correlation  rax(m)  as : 

rax(m) »  X^a(m+k)x(k) 

mcr  computes  :  Max{  raa(0)rxx(0))  } 

as  well  as  the  value  of  m  at  which  this  maximum  occurs. 

informations  required 

1 .1  Number  of  data  points  in  each  set  of  data. 

1 .2  Filenames,  starting  points  and  ending  points  of  the 
referenced  data  and  the  test  data. 

(vii)  fkusk 
description 

This  function  computes  the  Kurtosis  and  Skewness  ( 
see  (iv)  above  for  definitions )  of  a  power  spectrum.  The  input  to 
this  function  should  be  a  power  spectrum. 

informations  required 

Same  as  (iv)  except  that  x(t)  is  replaced  by  lX(f)l2  as  the  input. 


(viii)  sft 


description 

This  option  allows  the  user  to  view  all  the  computed  features 
on  the  screen  in  a  tabulated  form. 

informations  required 
None. 


4.2.6.  fsp  (fig  .4) 

description 

This  function  computes  the  power  spectrum  of  the  input  signal 
using  the  FFT  algorithm. 

informations  required 

1 .1  Filename  of  the  input  data. 

1 .2  Starting  point  and  ending  point  of  the  data. 

2.1  Number  of  spectral  points  required  ( power  of  2,  maximum 
is  512). 

2.2  FFT  order  { must  be  consistent  with  the  answer  to  2.1 ,  e.g., 
if  you  enter  256  in  2.1,  you  should  enter  8  here;  maximum 
is  9) 

3.1  An  output  filename  to  store  the  computed  magnitude 
spectrum. 

3.2  Number  of  output  points  to  be  stored. 


4.2.7.  hb  (fig.10) 

hb  stands  for  Hilbert  Transform;  the  real  part  and  imaginary  part 
of  an  analytic  signal  are  HiUrert  transforms  of  each  other. 

description 

This  function  computes  the  analytic  signal  corresponding  to 
the  input  using  signal  FFT. 

Let  s{t)  be  a  complex  signal  such  that 

S(t)-S/t)+/Si(t) 


1  2 


where  5^(1)  and  Si(t)  are  both  real  functions. 

If  s(t)  is  an  analytic  signal,  then  Si(t)  is  the  Hilbert  transform 
of  Sf(t)  or  s^(t)  is  the  negative  Hilbert  transform  of  Sj(t).  Another 
important  property  of  the  analytic  signal  s(t)  is  that  S(f)  =  0  for 
all  f<0. 

Our  implementation  of  this  function  is  based  on  this  latter 
property.  First  we  compute  X(k)  using  FFT.  Then  we  set  X(k)  *  0 
for  all  negative  frequencies  and  the  analytic  signal  is  then 
obtained  by  inverse  FFT.  The  magnitude  of  the  resulting  analytic 
signal  is  then  computed  and  stored.  This  signal  will  exhibit  the 
energy  distribution  of  the  input  in  the  time  domain. 

informations  required 

1 .1  Filename  of  the  input  signal. 

1 .2  Starting  point  and  ending  point  of  the  signal. 

2.1  Number  of  points  of  analytic  signal  required  (power  of  2, 
maximum  is  512). 

2.2  FFT  order  ( refer  to  (f)-2.2 ). 

3.1  An  output  filename  to  store  the  computed  analytic  signal. 

3.2  Number  of  output  points  need  to  be  stored. 


4.2.8.  nimp 

description 

Nonlinear  Mapping  is  an  algorithm  that  maps  vectors  from,  say, 
an  L-dimensional  space  to  some  lower  dimensional  space  ( for 
our  case,  this  is  usually  a  2  or  3  dimensional  space). 

Let  X,  be  vectors  in  the  originai  .nd  Y,  be  the 
corresponding  vectors  (i.e.  after  mapping  X|)  in  the  lower 
dimensional  space.  Let  denotes  the  Euclidean  distance 
between  X|  and  Xj  and  the  Euclidean  distance  between  Y| 
and  Yj.  This  algorithm  then  seeks  the  optimal  Y,'s  in  the  sense 
that 
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where 


C  -  Z,  df  , 

is  minimized. 

The  details  of  the  derivation  of  this  algorithm  is  described  in 
[3].  This  is  an  iterative  algorithm  such  that  iteration  stops 
when  E  is  less  than  some  prescribed  error  bound  supplied  by 
the  user. 

informations  required 

1 .1  Dimension  of  the  original  space  (  say  m  ). 

1 .2  Number  of  m-vectors  to  be  mapped. 

1 .3  Desired  dimension  of  the  space  to  be  mapped  onto. 

1 .4  A  filename  to  read  the  input  m-vectors  ( m  components  in 
each  input  vector ). 

2.1  Since  an  iterative  algorithm  is  used  in  this  part.  Iteration 
will  stop  when  either  the  number  of  iterations  reaches 
1000  or  the  error  measure  becomes  less  than  some  error 
bound  specified  by  the  user.  At  this  point,  the  user  will  be 
prompted  to  provide  an  errror  bound, 

3.1  After  a  successful  run,  the  user  will  be  prompted  to 
provide  an  output  filename  to  store  the  lower  dimensional 
vectors, 

4,2.9.  spi  (fig.13-16) 
description 

Spiking  filtering  [5]  is  a  deconvolution  operation  by  which  we 
estimate  the  transfer  function  of  a  defect  (which  we  model  as  a 
linear  time-invariant  system)  with  respect  to  some  reference. 

The  first  step  is  to  find  a  spiking  filter  h{k)  which  reduces  the 
reference  signal  (widely  known  as  a  wavelet  in  reflection 
seismology)  to  9(k-kg),  which  is  a  time-shifted  unit  sample,  i.e. 

h(k)  *  xjk)  -  9(k-k,) 

where  *  denotes  the  convolution  operator  and  kg  is  a  variable 
to  be  determined. 


Let  e  -  1/R  =  S.  (  h''(k).x„(k)  -  a(k-kj  f 

l\ 

where  h''(k)  is  a  least  mean  square  error  estimate  of  h(k).  Then  kg  is 
chosen  as  such  that  e  is  minimum  or  R  is  maximum.  R  is  known  as 
the  Performance  Index. 

Once  an  optimal  h''(k)  ( in  the  least  mean  square  error  sense ) 
has  been  obtained,  the  transfer  function  is  estimated,  up  to  a 
linear  time  shift  ko,  by  h''(k)*Xg(k). 


informations  required 

1 .1  Filename,  starting  point  and  ending  point  of 
signal. 

1 .2  Filename,  starting  point  and  ending  point  of 

1.3  A  filename  to  store  the  deconvolved  signal. 


the 

the 


reference 
test  signal. 


note:  spi  is  not  available  on  the  IBM-PC  version  of  this  package. 
Some  results  can  be  found  in  the  DEMO.exe  program. 


4.2.10.  wv  (fig.6) 

description 

Let  x(k)  be  the  input  signal,  0  <  k  <  N-1 .  The  wavelet 
transformation  X(/i,/3)  is  defined  as: 


X(Ai,5)  -  (1/V^)ZJ  x(k)cos((k-fl)/a)exp(-{(k-/^//i)2/2)} 

for  0^/?<N-1  : 

where  p  is  called  a  scale  factor. 


wv  will  compute  of  up  to  8  different  values  of  p 
and  display  all  transformed  wavelets  on  one  single  plot. 


From  our  experience,  although  the  pulse  echoes  from  hidden 
defects  of  different  geometries  may  themselves  look  very 
similar,  the  wavelet  transform  patterns  they  produced  are 


often  quite  different. 


informations  required 

1 .1  Filename,  starting  point  and  ending  point  of  the  input 
signal. 

1 .2  Number  of  data  points  in  the  input  wavelet. 

1 .3  Number  of  differnt  scales  required:  the  first  scale  and  the 
step  size  between  scales.  All  these  must  be  integers. 


4.2.11.  win 

descrip  tion 

This  is  another  deconvolution  operation  using  the  method  of 
Wiener  Filtering  [5]. 

The  Wiener  Filter  is  given  by 


H„(k)  -  X„*(k)/ {  |Xo(k)l2  +  1  /SNR(k)  } 


where  SNR(k)  =  Sx,(k)/ Snjk) 

is  the  ratio  of  the  power  spectral  density  of  Xp(k)  to  that  of  the 
noise. 

To  avoid  the  estimation  of  these  power  spectral  densities, 
SNR(k)  is  usually  replaced  by  a  small  positive  constant  1/q,  where 
q  is  called  the  noise  desensitizing  factor.  As  rule  of 
thumb,  q  is  taken  to  be  about  0.01  times  the  peak  value  of  |XQ(k)|2. 
The  required  transfer  function  is  then  estimated  by 

Inverse  FFT(  HJk)X(k) ) 


informations  required 

1 .1  Filename,  starting  point  and  ending  point  of  the  referenced 
signal. 

1 .2  Filename,  starting  point  and  ending  point  of  the  test  signal. 
2.1  Provide  a  noise  desensitizing  factor  q.  As  a  rule  of  thumb, 

q  can  be  chosen  to  be  about  0.01  times  the  maximum  of  the 
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magnitude  spectrum  of  the  referenced  signal. 

3.1  Number  of  spectral  points  to  be  used  in  the  filtering 
process  ( power  of  2,  maximum  is  512). 

3.2  FFT  order  (refer  to  (f)-2.2). 


4.2.12.  wind  (fig.11-12) 
description 

This  is  a  3-point  moving  average  operation  which  is  equivalent 
to  a  Hann's  frequency  window.  In  particular,  if  x(k)  is  the  input 
signal,  then  x(k)  will  be  replaced  at  the  output  by 

0.25{  x(k+1)  +  x(k-1) }  +  0.5x(k). 


required  informations 

1 .1  Filename,  starting  point  and  ending  point  of  the  input  data 
to  be  filtered. 

1 .2  Total  number  of  points  to  be  filtered. 

2.1  An  output  filename  to  store  the  filtered  signal. 

2.2  Number  of  output  points  need  to  be  stored. 


4.2.13.  fish  {fig.22) 
description 

Multiple  Fisher's  Discriminant  is  a  feature  space  compression 
algorithm.  The  details  for  Its  derivation  can  be  found  in  [6]. 

Suppose  at  the  beginning,  we  have  c  classes  (say  c,,C2.....  and 
each  of  P  patterns  in  an  n-dimensional  feature  space,  where  n  >  c. 
This  operation  can  optimally  map  these  n-vectors  into  ,for 
instance,  an  m-dimensional  space  such  that  m  <  c-1 . 

Let  X  be  one  of  the  available  patterns  of  known  class  in  the 
starting  feature  space.  We  need  to  find  an  optimal  linear 
transformation,  an  n  by  m  matrix  W  to  map  each  vector  x  into  a 
vector  y  in  the  lower  dimensional  space  such  that 


y  »  W^x 


First,  we  compute  the  within  class  and  between  class  scatter 
matrices  which  are  denoted  respectively  by  and  S^.  These 
matrices  are  defined  as  follows  : 

c 

Sv, » X  X  (x-m,){x-m,)T 

i=1  all  X  in  C; 

C 

Sb  =  X  P(m,-m){m,-m)T 

i=i 

where 

m,  =»  (1/P)  X  X 

all  X  in  Cj 

C 

m  =  f1/cPlX,Pm, 
ui 

The  matrix  W  is  found  such  that  the  criterion  function 

j(W)  -  IwTShWl/lW^s^wl 

is  maximized. 

The  solution  to  this  problem  can  be  obtained  by  solving 
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SbWj  - 


where  ^ij  is  known  as  the  generalized  eigenvalues  and  Wj  the 
corresponding  generalized  eigenvectors. 

It  can  be  easily  shown  that  the  rank  of  is  no  greater 
than  c-1  which  implies  that  no  more  than  c-1  of  the 
eigenvalues  are  non-zero.  Now,  let  m  be  the  number  of  non-zero 
eigenvalues  and  define 

W  =  {  Wg ...  } 

where  w^,  Wg, ...  are  the  eigenvectors  that  correspond  to 

those  non-zero  eigenvalues.  These  vectors  are  also  referred  to 
as  the  direction  vectors  for  obvious  reasons. 

Due  the  finite  precision  of  our  computers,  the  "unused" 
eigenvalues  will  not  be  exactly  zero.  Therefore,  in  practical 
implementation,  we  need  to  monitor  the  relative  sizes  of  all 
eigenvalues  and  determine  which  eigenvalues  should  be 
considered  as  zero.  ( From  our  experiences,  those  eigenvalues 
that  is  less  than  0.01  times  the  maximum  eigenvalue  can 
safely  be  taken  as  zero's.) 

informations  required 

1 .1  Number  of  clusters  of  training  patterns  (maximum  is  5). 

1 .2  Number  of  features  per  pattern  (maximum  is  6). 

1 .3  Number  of  training  patterns  per  cluster  (maximum  is  54). 

1 .4  A  filename  for  each  cluster  of  training  patterns. 

note:  We  have  assumed  that  each  cluster's  pattern  data  are 
stored  in  a  different  file.  For  example,  if  you  have  3 
clusters,  each  of  10  training  patterns  and  each  pattern 
consists  of  5  features,  you  need  to  store  these  data  in  3 
files,  and  each  file  contains  50  pieces  of  data. 

After  entering  the  above  informations,  the  following 
informations  will  be  dumped  on  the  screen: 

(i)  The  between-class  scatter  matrix  (S),) . 

(ii)  The  within-class  scatter  matrix  (S  J  . 

(iii)  The  generalized  eigenvalues  and  generalized 

eigenvectors  of  (S^,  S^). 

2.1  At  this  point,  the  user  is  prompted  to  provide  a  threshold 
to  choose  eigenvalues.  Usually,  this  threshold  is  less  than 
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0.01  times  the  maximum  eigenvalue.  Eigenvalues  less  than 
the  threshold  will  be  deleted  in  the  subsequent 
computations. 

3.1  If  no  numerical  ill-condition  is  encountered,  the  user  will 
be  prompted  to  provide  an  output  filename  to  store  the 
results  at  the  end  of  the  computations. 


4.2.14.  foly  (fig.21) 
description 

Foley  Sammon  Transformation  [8]  can  be  used  for  feature  space 
compression  when  there  are  only  two  classes  in  the  feature  space 
of  interest. 

Suppose  we  have  2  classes  of  patterns,  say  c^  and  Cg  in  an 
n-dimensional  feature  space  (the  starting  space)  and  we  want  to 
map  these  patterns  into  an  m-dimensional  space  (the  ending  sace) 
with  m  n.  Let  X  be  a  pattern  (or  vector)  in  the  starting  space  and 
y  its  correspondence  in  the  ending  space.  We  need  to  find  m 
n-vectors  ,  say  w,.  Wg . and  that 

y  *  W'^’x  where  W  =  { w, ,  Wg, ....  } 


These  vectors  are  called  direction  vectors. 

Define  the  criterion  function 

J(w)  -  (w'''(m,  -  mg))Vw''^S^w 

where  nrii  and  mg  are  the  mean  vectors  of  c,  and  Cg  respectively 
and  is  the  within  class  scatter  matrix  (see  4.2.13).  This  is  the 
the  same  criterion  function  used  in  Fisher's  linear  discriminant. 

The  first  direction  vector  is  computed  as  such  that  J(w^)  is 
maximized.  Since  our  primary  concern  is  the  direction  of  w^,  we 
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may  impose  an  additional  condition  that  »  1  .  The  solution 
is  found  as 

Wi  = 

6^2  =  dT(S^-i)2d 


where  d  =  (m^  -  m2) 


Then  Wg  is  found  also  by  maximizing  J(W2)  with  the  additional 
constraints  of  =1  and  Wg'^’w^^O  .  The  process  is  carried 
on  until  is  obtained  by  maximizing  J(w  J  with  additional 
constraints  that  »  1  and  that  =  0  for  all  i  ^  m. 

If  we  define  the  "discrim  values"  =  J(Wj),  we  shall  observe 
that 

ri>r2>....^r^ 


informations  required 

note:  The  user  should  first  provide  2  clusters  of  training 
patterns.  We  have  assumed  that  each  cluster  contains  the  same 
number  of  patterns  and  that  each  cluster’s  data  are  stored  in  a 
separate  file,  i.e.  2  input  files  are  required. 

1 .1  Number  of  training  patterns  per  cluster. 

1 .2  Number  of  features  per  pattern. 

1 .3  Desirable  number  of  compressed  features  per  pattern. 

1.4  A  filename  to  read  the  data  of  the  first  cluster. 

1 .5  A  filename  to  read  the  data  of  the  second  cluster. 

After  these  informations  are  entered,  the  following  results 
will  be  dumped  on  the  screen: 

(i)  The  within-class  scatter  matrix  Syy. 

(ii)  The  direction  vectors  Wj  . 

2.1  Number  of  test  patterns  ( whose  festures  are  to  be 
compressed ). 

2.2  A  filename  to  read  all  these  test  patterns. 

3.1  After  a  successful  computation,  the  user  will  be  prompted 
to  provide  a  filename  to  stored  the  compressed  features. 


21 


4.2.15.  kmean 


description 

Kmean  is  a  clustering  algorithm  which  iteratively  updates  the 
elements  of  each  cluster.  Suppose  there  are  P  test  patterns  needed 
to  be  grouped  into  K  clusters.  We  first  make  an  initial  guess  of  the 
K  cluster  centers  (each  cluster  center  is  the  mean  vector  of  all  the 
patterns  in  that  cluster).  If  no  a  priori  knowledge  is  available, 
these  initial  centers  may  be  chosen  arbitrarily  (for  instance,  the 
first  K  of  the  P  input  patterns).  Each  pattern  is  then  assigned  to  a 
cluster  such  that  the  Euclidean  distance  between  this  vector  and 
the  center  of  that  cluster  is  the  smallest.  When  all  patterns  have 
been  assigned,  each  cluster  center  is  updated  by  replacing  it  with 
the  mean  vector  of  the  cluster.  Iteration  goes  on  until  all  cluster 
centers  stablize. 


informations  required 

1 .1  Number  of  test  patterns  to  be  classified. 

1 .2  Number  of  features  in  each  pattern. 

1 .3  A  filename  to  read  all  these  test  patterns. 

2.1  Desired  number  of  clusters  (K). 

2.2  Guess  the  positions  of  the  initial  centers  of  each  cluster  in 
the  input  file. 


After  these  informations  are  entered,  conclusions  will  be 
dumped  on  the  screen. 

note:  These  results  will  not  be  stored  in  a  file.  The  user  may 
want  to  dump  these  on-screen  results  onto  a  printer  using  a 
DOS  command. 
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4.3  Wig.exe 

WIG.exe  is  another  executable  file  that  computes  the  Wigner-Ville 
Distribution  [11]  of  an  input  signal  and  will  provide  a  3-D  graphics 
display  of  the  results.  Two  graphics  at  different  view  point  will  be 
displayed  (one  after  the  other). 


description 

Let  x(t)  be  the  input  signal,  the  Wigner-Ville  distribution  is 
defined  as; 

-  J  x{t+/3/2)  x‘(t-;?/2)  exp(-jnt) 

where  x*(t)  denotes  the  complex  conjugate  of  x(t). 

This  is  a  time-frequency  distribution  and  is  believed  to  be  capable 
of  providing  time-domain  and  frequency-domain  informations 
simultaneously.  Such  a  property  is  very  useful  for  defect  classfication 
problems  since  both  time-domain  and  frequency-domain  both  provide 
relevant  informations . 


informations  required 

1 .  A  filename  to  read  the  input  signal. 

2.  Starting  point  (ist)  and  Ending  point  (ien)  of  this  signal. 

note:  Only  a  maximum  of  64  points  are  allowed  in  the  input 
data.  Hence,  we  require  that :  (ien  -  ist)  <  64  . 

3.  Number  of  time  points  for  which  Wigner-Ville  distribution  is  to 
be  computed  (  maximum  is  64  ). 
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5.  Computer  Results 


In  this  section,  we  present  some  results  that  are  obtained  through  the 
use  of  this  package. 

The  reference  signal  r(t)  and  the  test  signal  x(t)  shown  in  figs.(1) 
and  (2)  are  pulse  echoes  obtained  by  using  a  15  MHz  ultrasonic  transducer 
and  digitized  with  a  sampling  rate  of  100  MHz.  Both  signals  have  a  record 
length  of  512.  r(t)  is  the  pulse  echo  from  a  flawless  aluminium  block 
while  x(t)  is  the  pulse  echo  from  a  block  with  a  hidden  angular  cut. 

Fig. (3)  shows  the  Power  Spectrum  of  r(t)  computed  by  a  512-point 
FFT  (Fast  Fourier  Transform).  The  frequency  resolution  is  512  points  per 
100  MHz. 

Fig.(4)  is  the  cross-correlation  function  Rpj(fi)  ^or  *  -1500  ns  to 
1500  ns  (note:  this  function  has  been  shifted  by  1500  ns  so  that  the  time 
origin  is  at  zero). 

Fig. (5)  shows  the  (AR)  Power  Spectrum  of  r(t)  computed  by  Burg's 
Technique.  The  AR  order  is  50  and  the  frequency  resolution  is  512  points 
per  100  MHz. 

The  signal  h^(t)  shown  is  fig. (6)  is  the  impulse  response  extracted 

from  x(t)  by  Wiener  filtering  with  r(t)  as  the  reference  signal.  The  noise 
desensitizing  factor  used  is  equal  to  0.01  times  the  peak  power  shown  in 
fig.(3). 


Fig. (7)  is  the  result  of  time-averaging  h^(t). 

Fig. (8)  shows  the  analytic  signal  of  hj^{t). 

In  fig. (9),  we  have  simulated  a  pulse  echo  from  a  2  layer  composite 
material  by  concatenating  2  single  reflections.  We  call  this  signal  the 
composite  signal  c(t).  The  impulse  response  hj,(t)  extracted  from  this 
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signal  by  Wiener  Filtering  { with  r(t)  as  reference )  is  shown  in  fig.{10). 
Fig.(1 1)  is  a  time-averaged  version  of  hg(t)  while  fig. (12)  is  the  analytic 
signal  of  h^{X). 

In  figs.(13)  to  (17)  we  demonstrate  the  different  steps  of  a  spiking 
filtering  operation  (note:  the  intermediate  results  shown  here  are  not  all 
accessible  by  the  user  of  this  package).  First,  let  the  reference  signal  (or 
the  reference  wavelet)  be  ref(t)  as  shown  in  fig. (13).  The  spiking  filter 
(or  inverse  filter)  h3(t)  shown  in  fig.(14)  is  computed  such  that  the 

performance  index  is  maximum  in  the  performance  index  profile  shown  in 
fig. (15).  We  may  check  the  accuracy  of  this  spiking  filter  by  convolving 
ref(t)  with  hg(t)  and  the  result  of  this  convolution  is  shown  in  fig.(16). 

Now,  let  the  test  signal  of  interest  be  s(t)  which  is  made  up  of  2 
overlapped  wavelets,  such  that  s(t)  «  ref(t)  +  0.5ref(t-5).  The  final 
result,  which  is  given  by  s(t)*hg(t),  is  shown  in  fig.(17).  You  may  observe 

that  the  occurence  of  the  two  spikes  in  fig.(17)  corresponds  exactly  to 
the  occurence  of  the  two  wavelets  in  s(t).  Fig. (18)  is  another  example  of 
spiking  filtering  in  which  the  test  signal  is  T15A1  and  the  reference 
signal  is  T15A0  ;  these  waveforms  can  be  found  in  Appendix  B  of  this 
report. 

Figs.(19)  and  (20)  show  the  Wigner-Ville  distribution  of  x(t).  The 
numbers  of  frequency  points  and  time  points  are  both  equal  to  64. 

Fig.(21)  and  (22)  are  scatter  plots  which  show  the  results  of 
feature  space  compression  from  a  4-dimensional  to  a  2-dimensional 
space  by  Foley-Sammon  Transformation  and  Multiple  Fisher  Discriminant 
respectively.  In  these  figures,  f1  and  f2  denote  the  2  dimensions  after 
the  compression.  The  original  data  used  are  Iris  data  which  are  tabulated 
in  Table  1  on  page  32  of  this  report.  The  2  classes  of  patterns  (Iris)  we 
chose  for  this  demonstration  are  "versicolor"  and  "setosa";  50  patterns 
from  each  class  have  been  used. 


power  power  ^  amplitude 

too  800  (S'  1.0*+002  t.0«-f007  8.0a-f007  il:  -096  0  600  1800 


Reference:r(t) 


time(zlO^s) 

.  1  A  reference  signal  r(t) 


FFT_Spectrum 


frequency(xO.  195_MKz) 

.  3  Power  spectrum  of  r(t)  by  FFT 


Burg’s_Spectrum 

r  -r  r  r  “i — »  -r  ~  i  ’»"n — r-t—r—r— — r-T — r 


frequency(xO.  195_lIHz) 


Fig.  5  Power  spectrum  of  r(t)  by  Burg's 
technique. 


time(zlO_na) 

Fig.  2  A  test  signal  x(t) 


ilme(zl0_n8) 


Fig.  4  Correlation  function 


Wiener_Filter(x/r) 

—  T~  r-i  -t  T  -1 — r— r"*! — 9-  T  i—r-i  -1  y 


tlme(zl0_n8) 

Fig.  6  hjj(t),  impulse  response  by 
Wiener  filtering. 
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Smoothing  (x/r) 


Anal3’-tic_Signal(x/ r) 


Fiy.  8  Analytic  signal  of  ^^(t) 


Composite;c(t) 


Smoothing(c/r) 


'\Viener_Filter(c/r) 


tlme(xl0_n8) 


Fig.  10  h^(t),  impulse  response  by  Wiener 

filtering 
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Composite— Signal  Performance— Index  reference 


Fi^'ures  13-17  show  steps  of  spihia^  filtering. 


Spiking-Filter 


Spiking -Filter 
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Figure  21  T\yO  clusters  of  Iris  data  set  obtained  by 
Foley -Saitiraon  algorithm 
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Multiple— Fisher  (irse,ve) 


fi 


Figure  22  Two  clusters  of  Iris  data  set  obtained 

by  multiple  Fisher  discriminant  algorithm. 
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Table  1 


List  of  Iris  Data  Set 


iris  setosa 

/m  versicolor 

iris  virginica 

Sepal 

length 

Sepal 

width 

Petal 

length 

Petal 

width 

Sepal 

length 

Sepal 

width 

Petal 

length 

Petal 

width 

Sepal 

length 

Sepal 

width 

Petal 

length 

Petal 

width 

5.1 

3.5 

1.4 

0.2 

7.0 

3.2 

4.7 

1.4 

6.3 

3.3 

6.0 

2.5 

4.9 

3.0 

1.4 

0.2 

6.4 

3.2 

4.5 

1.5 

5.8 

2.7 

5.1 

1.9 

4.7 

3.2 

1.3 

0.2 

6.9 

3.1 

4.9 

1.5 

7.1 

3.0 

5.9 

2.1 

4.6 

3.1 

1.5 

0.2 

5.5 

2.3 

4.0 

1.3 

6.3 

2.9 

5.6 

1.8 

5.0 

3.6 

1.4 

0.2 

6.5 

2.8 

4.6 

1.5 

6.5 

3.0 

5.8 

2.2 

5.4 
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6.  Concluding  Remarks 
and  Future  Developments 


The  lUNDE  package  is  developed  with  the  primary  intentions  of 
being  very  easy  to  use  and  that  the  user  can  familiarize  himself  with  all 
these  operations  very  quickly.  It  is  our  belief  that  this  package  can  serve 
as  a  very  effective  tool  that  assists  engineers  and  researchers  in 
practical  applications  to  NDE  problems. 

The  IBM-PC  (or  compatibles)  version  of  this  package  has  a  very 
powerful  graphics  capability.  Many  different  kinds  of  peripherals  are 
allowed.  The  user  can  run  the  EQUIP.exe  program  to  see  the  details. 

A  second  version  of  the  lUNDE  package  is  being  developed  by 
Information  Research  Lab.,  Inc.  and  is  aimed  towards  the  following  goals  ; 

1.  Expansion  of  the  signal  processing  library  by  including  such 
functions  as  Spectral  Zooming  and  Spectral  Expansion, 
Bicorrelational  Bispectral  Analysis  and  Model-Based  Bispectral 
Analysis,  Spectral  Extrapolation  and  L2  Deconvolutions,  etc. 

2.  Integrate  with  the  RCE™  Neual  Network  developed  by  Nestor®  Inc. 
to  perform  more  reliable  and  versatile  pattern  classifications. 

3.  Include  the  ability  to  "talk"  with  some  specific  ultrasonic 
hardwares  such  as  pulser/receiver  and  high  speed  digitizer  to  build 
up  an  on-line  defect  classification  system. 

4.  Expansion  of  the  graphic  display  fuctions  by  adding  the  options  for 
3-dimensional  surface  plots,  contour  plots  and  multi-window 
graphic  displays. 

5.  Enhancement  of  the  user  interface:  the  new  version  lUNDE  is  all 
menu-driven  with  pop-up  dialog  box  which  adds  to  the  use  of  this 
package  a  lot  more  fun. 
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Appendix  B  :  Ultrasonic  Waveforms 


We  have  included  some  ultrasonic  pulse  echo  samples  in  this  package.  All 
these  signals  are  sampled  with  a  sampling  rate  of  100  MHz.  The  names  of 
these  signals  should  be  interpreted  as  follows  : 

(1) T5  or  T15  indicates  the  type  of  ultrasonic  transducers  used  for 
acquiring  these  signals.  T5  indicates  a  5  MHz  transducer  and  T1 5 
indicates  a  15  MHz  transducer. 

(2) A,  B.  C  and  D  indicate  the  relative  sizes  of  the  hidden  defects  in 
the  aluminium  blocks  from  which  these  pulse  echoes  are  obtained; 
such  that 

sizes  of  A  :  B  :  C:  D  =  8  : 4  : 2  : 1 

(3)  AO  :  no  hidden  defects  ( a  flawless  sample ). 

A1 ,  B1 ,  Cl  and  D1  :  hidden  defects  are  flat-topped  cuts. 

A2,  B2,  B3,  C2,  C3. 

D2,  D3  and  D4  :  hidden  defects  are  angular  cuts. 

A3,  B4,  C4  and  D5  :  hidden  defects  are  circular  holes. 
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